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ABSTRACT 

Aims. We revisit the study of the mass functions and the bias of dark matter halos. 

Methods. Focusing on the limit of rare massive halos, we point out that exact analytical results can be obtained for 
the large-mass tail of the halo mass function. This is most easily seen from a steepest-descent approach, that becomes 
asymptotically exact for rare events. We also revisit the traditional derivation of the bias of massive halos, associated 
with overdense regions in the primordial density field. 

Results. We check that the theoretical large-mass cutoff agrees with the mass functions measured in numerical simula- 
tions. For halos defined by a nonlinear threshold S — 200 this corresponds to using a linear threshold Sl — 1.59 instead 
of the traditional value ~ 1.686. We also provide a fitting formula that matches simulations over all mass scales and 
obeys the exact large-mass tail. Next, paying attention to the Lagrangian-Eulerian mapping (i.e. corrections associated 
with the motions of halos), we improve the standard analytical formula for the bias of massive halos. We check that 
our prediction, which contains no free parameter, agrees reasonably well with numerical simulations. In particular, it 
recovers the steepening of the dependence on scale of the bias that is observed at higher redshifts, which published 
fitting formulae did not capture. This behavior mostly arises from nonlinear biasing. 

Key words. Cosmology: large-scale structure of Universe; gravitation; Methods: analytical 



1. Introduction 

The distribution of nonlinear virialized objects, such as 
galaxies or clusters of galaxies, is a fundamental test of cos- 
mological models. First, this allows us to check the validity 
of the standard cosmological scenario for the formation of 
large-scale structures, where nonlinear objects form thanks 
to the amplification by gravitational instability of small pri- 
mordial density fluctuations, built for instance by an early 
inflationary stage (e.g., Peebles 1993; Peacock 1998). For 
cold dark matter (CDM) scenarios (Peebles 1982), where 
the amplitude of the initial perturbations grows at smaller 
scales, this gives rise to a hierarchical process, where in- 
creasingly large and massive objects form as time goes on, 
as increasingly large scales turn nonlinear. This process has 
been largely confirmed by observations, which find smaller 
galaxies at very high redshifts (e.g., Trujillo et al. 2007) 
while massive clusters of galaxies (which are the largest 
bound objects in the Universe) appear at low redshifts 
(e.g., Borgani et al. 2001). Second, on a more quantitative 
level, statistical properties, such as the mass function and 
the two-point correlation of these objects, provide strong 
constraints on the cosmological parameters (e.g. through 
the linear growth factor D + (t) of density perturbations) 
and on the primordial fluctuations (e.g. through the ini- 
tial density power spectrum Pj,(fc)). For these purposes, 
the most reliable constraints come from observations of the 
most massive objects (rare-event tails) at the largest scales. 
Indeed, in this regime the formation of large-scale struc- 
tures is dominated by the gravitational dynamics (bary- 
onic physics, which involves intricate processes associated 
with pressure effects, cooling and heating, mostly occurs at 
galactic scales and below), which further simplifies as one 
probes quasi-linear scales or rare events where effects asso- 



ciated with multiple mergings can be neglected. Moreover, 
in this regime astrophysical objects, such as galaxies or clus- 
ters of galaxies, can be directly related to dark matter halos, 
and their abundance is highly sensitive to cosmological pa- 
rameters thanks to the steep decline of the high-mass tail 
of the mass function (e.g., Evrard 1989). 

Thus, the computation of the halo mass function (and 
especially its large-mass tail) has been the focus of many 
works, as it is one of the main properties measured in galaxy 
and cluster surveys that can be compared with theoretical 
predictions. Most analytical derivations follow the Press- 
Schechter approach (Press & Schechter 1974; Blanchard et 
al. 1992) or its main extension, the excursion set theory of 
Bond et al. (1991). In this framework, one attempts to esti- 
mate the number of virialized objects of mass M from the 
probability to have a linear density contrast Sl at scale M 
above some given threshold 5 C . Thus, one identifies current 
nonlinear halos from positive density fluctuations in the 
initial (linear) density field, on a one-to-one basis. This is 
rather well justified for rare massive objects, where one can 
expect such a link to be valid since such halos should have 
remained well-defined objects until now (as they should 
have suffered only minor mergers). By contrast, small and 
typical objects have experienced many mergers and should 
be sensitive to highly non-local effects (e.g. tidal forces, 
mergers), so that such a direct link should no longer hold, 
as can be checked in numerical simulations (Bond et al. 
1991). As noticed by Press & Schechter (1974), the simplest 
procedure only yields half the mass of the Universe in such 
objects (essentially because only half of the Gaussian initial 
fluctuations have a positive density contrast, whatever the 
smoothing scale), which they corrected by an ad- hoc mul- 
tiplicative factor 2. In this respect the main result of the 
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excursion set theory was to provide an analytical derivation 
of this missing factor 2, in the simplified case of a top-hat fil- 
ter in Fourier space (Bond et al. 1991). Then, it arises from 
the fact that objects of mass larger than M are associated 
with configurations such that the linear density contrast 
goes above the threshold S c at some scale M' > M, which 
includes cases missed by the Press-Schechter prescription 
where the linear density contrast decreases below this scale 
M' so that 6l < S c at scale M ( "cloud-in-cloud" problem). 
The characteristic threshold 5 C is usually taken as the lin- 
ear density contrast reached when the spherical collapse dy- 
namics predicts collapse to a zero radius. In an Einstein-dc 
Sitter universe this corresponds to S c ~ 1.686 and to a non- 
linear density contrast 6 ~ 177 (assuming full virialization 
in half the turn- around radius). This linear threshold only 
shows a very weak dependence on cosmological parameters. 

Numerical simulations have shown that the Press- 
Schechter mass function (PS) is reasonably accurate, es- 
pecially in view of its simplicity. Thus, it correctly pre- 
dicts the typical mass scale of virialized halos at any red- 
shift, as could be expected since the large-mass behavior 
is rather well justified. In addition, it predicts a universal 
scaling that appears to be satisfied by the mass functions 
measured in simulations, that is, the dependence on halo 
mass, redshift and cosmology is fully contained in the ratio 
v = S c /a(M, z), where a(M,z) is the rms linear density 
fluctuation at scale M. However, as compared with numer- 
ical results it overestimates the low-mass tail and it under- 
estimates the high-mass tail. This has led to many numer- 
ical studies which have provided various fitting formulae 
for the mass function of virialized halos, written in terms 
of the scaling variable v (or a) (Sheth & Tormcn 1999; 
Jenkins et al. 2001; Reed et al. 2003; Warren et al. 2006; 
Tinker et al. 2008). We may note here that a theoretical 
model that attempts to improve over the PS mass func- 
tion is to consider the ellipsoidal collapse dynamics within 
the excursion-set approach, to take into account the devia- 
tions from spherical symmetry for intermediate mass halos 
(Sheth et al. 2001). Note that, as described in the original 
paper (and emphasized in Robertson et al. 2009), at large 
mass this would recover the spherical collapse. However, the 
halo mass obtained by such methods is generally underes- 
timated for non center-of-mass particles (since analytical 
computations assume that particles are located at the cen- 
ter of their halo, i.e. they only consider the linear densities 
within spherical cells centered on the point of interest). To 
correct for this effect, in practice one treats the threshold 
5 C as a free parameter, close to 1.6, to build fitting formulae 
from numerical simulations. For instance, this correction is 
contained in the parameter a in the exponential cutoff of 
Sheth et al. (2001). 

A second property of virialized halos that can be used 
to constrain cosmological models, beyond their number 
density, is their two-point correlation. Indeed, observations 
show that galaxies and clusters do not obey a Poisson dis- 
tribution but show significant large-scale correlations (e.g., 
McCracken et al. 2008; Padilla et al. 2004). In particular, 
their two-point correlation roughly follows the underlying 
matter correlation, up to a multiplicative factor b 2 , called 
the bias, that grows for more massive and extreme objects. 
Following the spirit of the Press-Schechter picture, Kaiser 
(1984) found that this behavior arises in a natural fashion if 
halos are associated with large overdensities in the Gaussian 
initial (linear) density field, above the threshold S c . This 



was further expanded by Bardeen et al. (1986) and Bond & 
Myers (1996), who considered the clustering of peaks in the 
Gaussian linear density field. A simpler derivation, based on 
a peak-background split argument, and taking care of the 
mapping from Lagrangian to Eulerian space, was presented 
in Mo & White (1996). It provides a prediction for the bias 
b(M) as a function of halo mass, in the limit of large dis- 
tance, r — > oo, that agrees reasonably well with numerical 
simulations. However, as for the PS mass function, in order 
to improve the agreement with numerical results various fit- 
ting formulae have been proposed (Sheth & Tormen 1999; 
Hamana et al. 2001; Pillepich et al. 2009). Again, since the 
ellipsoidal collapse model reduces to the spherical dynam- 
ics for rare massive halos, it also requires free parameters 
to improve its accuracy, but the latter are consistent with 
those used for the mass function (Sheth et al. 2001). 

In this article we revisit the derivation of the mass 
function and the bias of rare massive halos, following the 
spirit of the Press & Schechter (1974) and Kaiser (1984) 
approaches. That is, we use the fact that large halos can 
be identified from overdensities in the Gaussian initial (lin- 
ear) density field. First, we briefly review in section [5] some 
properties of the growth of linear fluctuations and of the 
spherical dynamics in ACDM cosmologies. Next, we recall 
in section [3] that in the quasi-linear regime (i.e. at large 
scales), the probability distribution of the nonlinear density 
contrast 5 r within spherical cells of radius r can be obtained 
from spherical saddle-points of a specific action <S, for mod- 
erate values of S r where shell-crossing does not come into 
play. We also discuss the properties of these saddle-points as 
a function of mass, scale and redshift. Then, we point out in 
section[l]that this provides the exact exponential tail of the 
halo mass function. This applies to any nonlinear density 
contrast threshold S that is used to define halos, provided 
it is below the upper bound <5 + where shell-crossing comes 
into play. We compare our results with numerical simula- 
tions and we give a fitting formula that applies over all mass 
scales and satisfies the exact large-mass cutoff. Next, we re- 
call in section [5] that these results also provide the density 
profile of dark matter halos at outer radii (i.e. beyond the 
virial radius) in the limit of large mass. Finally, we study 
the bias of massive halos in section [6l paying attention to 
some details such as the Lagrangian-Eulerian mapping, and 
we compare our results with numerical simulations. We con- 
clude in section [7l 



2. Linear perturbations and spherical dynamics 

We consider in this article a flat CDM cosmology with 
two components, (i) a non-relativistic component (dark 
and baryonic matter, which we do not distinguish here) 
that clusters through gravitational instability, and (ii) an 
uniform dark energy component that does not cluster at 
the scales of interest, with an equation-of-state parameter 
w = Pdc/pdc- For the numerical computations we shall focus 
on a ACDM cosmology, where the dark energy is associated 
with a cosmological constant that is exactly uniform with 
w = — 1. However, our results directly extend to curved 
universes (i.e. fik ^ 0) and to dark energy models with a 
possibly time- varying w(z), as long as we can neglect the 
dark energy fluctuations on the scales of interest, which is 
valid for realistic cases. Focussing on the case of constant w, 
we first recall in this section the equations that describe the 
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dynamics of the background and of linear matter density 
perturbations, as well as the nonlinear spherical dynamics. 

The evolution of the scale factor a(t) is determined by 
the Friedmann equation (Wang & Steinhardt 1998), 



H 2 (t) 
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with H(t) = -, (1) 



where subscripts denote current values at z = 0, when 
a = 1, and a dot denotes the derivative with respect to 
cosmic time t. On the other hand, the density parameters 
vary with time as 
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Next, introducing the matter density contrast, <5(x, t) = 
(p m (x, t) — p m )/~p m , where /5 m (t) is the mean matter density, 
linear density fluctuations grow as 



mass M in a uniform universe with the same cosmology. 
This also implies that the density p m writes as 

-3 



Pm(t) = p m (t)y(t)~ 
Then, Equation J5|) leads to 
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Of course, we can check that in the linear regime, where 
Ul = 1 — Sl/3, we recover Eq.Q for the linear growth 
of 5l- Then, to obtain the nonlinear density contrast, 
8{z) — p m /~p m — 1, associated with the linear density con- 
trast, 5l(z), at a given redshift z, we solve Eg. ([Til) with 
the initial condition y(zi) = 1 — 8^/3 and y'{z{) — —Sli/3 
at some high redshift Z[, with 8lJ8l — D + (zi) / D + (z). For 
any redshift z this defines a mapping <5l M> <5 = J-(5l) that 
fully describes the spherical dynamics before shell-crossing. 



2 Q -k 

a 



Angp m 5 L = 0, 



(3) 



where the subscript L denotes linear quantities. For nu- 
merical purposes it is convenient to use the logarithm of the 
scale factor as the time variable. Then, using the Friedmann 
equation (TTJ) the linear growth factor D + {t) evolves from 
Eq.© as 



D'l 



1 3 o 



D' 



-n m D+ = o, 



(4) 



where we note with a prime the derivative with respect 
to In a, as D' + = d-D+/dlna. Then, the normalized linear 
growth factor, g(t), defined as 

5( t ) = ^77r and s(* = o) = i, 



a(t) 
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5 3 o 



.9 



-(1 - w)fl dc g = 0, 
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with the initial conditions g — > 1 and g' — > at a — > 0. 

In the following we shall also need the dynamics of 
spherical density fluctuations. For such spherically symmet- 
ric initial conditions, the physical radius r(i), that contains 
the constant mass M until shell-crossing, evolves as 



4tt£ 



r [p m + (1 + 3w)p dc ] with p n 
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4 7rr 3 



(7) 



Note that p m is the mean density within the sphere of radius 
r (and not the local density at radius r). Using again the 
Friedmann equation (TTJ) this reads as 



p m + (1 + 3w)p dc 

2 (Pm + Pde) 



r = 0. 



(8) 



As for the linear growth factor D + (t), it is convenient to 
introduce the normalized radius y(t) defined as 



y(t) = 



?(*) 



and q(t) oc a(t), y(t = 0) = 1. 



(9) 



Thus, q(t) is the Lagrangian coordinate of the shell r(t), 
that is, the physical radius that would enclose the same 




Fig. 1. The function ^(Sl) that describes the spherical dy- 
namics when there is no shell-crossing, at z = 0. The solid 
line corresponds to the ACDM cosmology (with J7 m o = 
0.27) and the dashed line to the Einstein-de Sitter case 
fi m o = 1. The second peak corresponds to a second collapse 
to the center, but for our purposes we only need T(bi) be- 
fore first collapse. 

We compare in Fig. [JJ the function F{8l) obtained at 
z = within the ACDM cosmology that we consider in this 
article (solid line) with the Einstein-de Sitter case (dashed 
line), where it has a well-known parametric form (Peebles 
1980). As is well known, we can check that the dependence 
on f2 m o is very weak until full collapse to a point, which 
occurs at slightly lower values of 8l for low O m o. We show 
in Fig. [2] the linear density contrasts, 8l = J r_1 (<5), asso- 
ciated with three nonlinear density contrasts, 8 = 100, 200 
and 300, as a function of the cosmological parameter f2 m (z). 
In agreement with Fig. [TJ they show a slight decrease for 
smaller f2 m . For 8 — 200 a simple fit (dashed line) is pro- 
vided by 

5 = 200: <5 L ~ 1.567 + 0.032 (f2 m + 0.0005) ' 24 , (12) 

which agrees with the exact curve to better than 5 x 10~ 4 
for n m > 0.2, which covers the range of practical interest. 

In this article we consider Gaussian initial fluctuations, 
which arc fully defined by the linear density power spectrum 



(5 L (k 1 )^(k 2 )) = 6 D {ki + k 2 ) P L (h), 



(13) 
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Fig. 2. The linear density contrast 6l = J 7 ^ 1 ^) associated 
with the nonlinear density contrast 8, through the spheri- 
cal dynamics, as a function of the cosmological parameter 
Q m (z) at the redshift of interest. We show the three cases 
8 = 100, 200 and 300. The dashed line is the fit jT2j to the 
case 5 = 200. 

where 8d is the Dirac distribution and we normalized the 
Fourier transform as 



(S(x) = / dke ikx <5(k), 



(14) 



where x and k are the comoving spatial coordinate and 
wavenumber. This gives for the linear density two-point 
correlation 



C L (xi,x 2 ) = (S L (xi)5 L (x 2 )) 



&|x 2 -xi| 



(15) 



We also introduce the smoothed density contrast, 5 r (x), 
within the sphere of radius r and volume V around position 
x, 



5 r (x) 



— <5(x + x') = / dke ik x S(k) W{kr), (16) 
v V 



with a top-hat window that reads in Fourier space as 



W{kr) = I — e ikx = 3 



sin(fcr) — fcrcos(fcr) 



(17) 



Then, in the linear regime, the cross-correlation of the 
smoothed linear density contrast at scales r\ and r 2 and 
positions xi and x 2 = xi + x reads as 

O'n.raC^) = (^i>i(xi)5ir 2 (xi + x)) 

= An J dk k 2 P L (k)W(k ri )W(kr 2 f m ^ X \ (18) 

In particular, ay — oy jr (0) is the usual rms linear density 
contrast at scale r. 



3. Distribution of the density contrast 

We recall here that in the quasi-linear limit, ay — ► 0, the 
distribution V(8 r ) of the nonlinear density contrast 8 r at 



scale r can be derived from a steepest-descent method 
(Valageas 2002a). In agreement with the alternative ap- 
proach of Bernardeau (1994a), this shows that rare-event 
tails are dominated by spherical saddle-points, which we use 
in sections 0HS1 to obtain the properties of massive halos. 

Since the system is statistically homogeneous we can 
consider the sphere of radius r centered on the origin x = 0. 
Then, we first introduce the cumulant generating function 

<p(y), 



e - V (y)/a r = 



yS r /cr T 



(19) 



which determines the distribution V[8 T 
verse Laplace transform 
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through the in- 



(20) 



In Ea. lflUl) we rescaled the cumulant generating function by 
a factor a% so that it has a finite limit in the quasi-linear 
limit ay — > for the Gaussian initial density fluctuations 
that we study in this article (Bernardeau et al. 2002). In 
particular, its expansion at y — reads as 



<p(y) 
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p=2 



i-y) p (%)c 



pl a 2 r {p - 1] ' 



(21) 



Then, the average (IT91 can be written as the path integral 



,-<p(v)/<r r 



■DS L (q) e 



-S[S L ]/o 



(22) 



where C L 1 is the inverse matrix of the two-point correlation 
(ITSI) and the action S reads as 



S\8r 



y5r[5L] + ^5 L .Cl\5 L 



(23) 
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Here <5,[<5l] is the nonlinear functional that affects 
to the initial condition defined by the linear density 
field <5l(q) the nonlinear density contrast S r within 
the sphere of radius r, and we note Sl-C^ .8l = 
/ dqidq 2 (5i(qi)C i ^ 1 (qi, q 2 )5i(q 2 ). Note that the action S 
does not depend on the normalization of the linear power 
spectrum since both of and Cl are proportional to Pl- 

In the quasi-linear limit, a r — > 0, as shown in Valageas 
(2002a) the path integral (|2"2"|) is dominated by the mini- 
murrQ of the action S[Sl], 



: <p(y) -t min S[S L ] 
Mq) 



(24) 



Using the spherical symmetry of the top-hat window W, in 
agreement with the approach of Bernardeau (1994a), one 
obtains a spherical saddle-point with the radial profile 



(25) 



where a qtg / = cr q<q i{Q, 0) and q is the Lagrangian coordinate 
associated with the Eulerian radius r through the spher- 
ical dynamics recalled in section [21 and q' is a dummy 



1 As described in details in Valageas (2002a), depending on 
the slope of the linear power spectrum the saddle-point ()25|l 
may only be a local minimum or maximum, but it still governs 
the rare-event limit, in agreement with physical expectations. 
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O 2 4 6 8 lO 



q'/q 

Fig. 3. The radial profile ([2"5jl of the linear density contrast 
dijqi of the saddle point of the action 5[<5l(q)]- We show the 
profiles obtained with a ACDM cosmology for the masses 
M = 10 n ,10 12 ,10 13 ,10 14 and lO^/i^Mg. A larger mass 
corresponds to a lower ratio Sl^ l&Lq at large radii q'/q > 1. 

Lagrangian coordinate along the radial profile (hereafter 
we denote by the letter q Lagrangian radii, which are as- 
sociated with the linear density field, whereas r denotes 
Eulcrian radii associated with the nonlinear density con- 
trast 5 r ). Thus, the radii q' and r 1 are related by 

q' 3 = (1 + S r >) r' 3 , with 5 r > = J~(5 ) , (26) 

where the function S r i = F{5Lq'), that describes the spher- 
ical dynamics, was obtained below Eq. (|ll[) and shown in 
Fig.ffl 

We can note that the profile ([25]) is also the mean con- 
ditional profile of the linear density contast SLq', under the 
constraint that is it equal to a given value SLq at a given 
radius q (e.g., Bernardeau 1994a). The reason why the non- 
linear dynamics gives back this result is that the nonlinear 
density contrast S r only depends on the linear density con- 
trast S Lq at the associated Lagrangian radius q, through the 
mapping S r = F{SL q )- Indeed, as long as shell-crossing does 
not modify the mass enclosed within the shell of Lagrangian 
coordinate q, its dynamics is independent of the motion of 
inner and outer shells (thanks to Gauss theorem) . Then, in 
order to obtain the minimum of the action S we could pro- 
ceed in two steps. First, for arbitrary Lagrangian radius q 
and linear contrast SLq, we minimize S with respect to the 
profile SLq' over q' ^ q. From the previous argument, only 
the second Gaussian term in ([23)) varies so that this partial 
minimization leads to the profile (1231) (indeed, for Gaussian 
integrals the saddle-point method is exact). Second, we 
minimize over the Lagrangian radius q (or equivalently over 
Sl or S r ), which leads to Eq.([2T)|) below. Here we also use the 
fact that a spherical saddle-point with respect to spherical 
fluctuations is automatically a saddle-point with respect to 
arbitrary non-spherical perturbations and it can be seen 
that for small y one obtains a minimum (then we assume 
that at finite y strong deviations from spherical symmetry 
do not give rise to deeper minima, which seems natural from 
physical expectations). We refer to Valageas (2002a, 2009b) 
for more detailed derivations. 

Note that the shape of the linear profile (|2"5j) depends 
on the shape of the linear power spectrum, whence on the 



mass scale M of the saddle point for a curved CDM linear 
power spectrum, but not on redshift. We show in Fig.[3]the 
profile ([2"5")l obtained for several masses M. For a power-law 
linear power spectrum, of slope n, Ea. (|2"5f leads to SLq' oc 
q/-(n+3) a ^ i ar g e radii^ q' — > qq. Then, since for a CDM 
cosmology n increases at larger scales, the profile shows a 
steeper falloff at large radii for larger mass, in agreement 
with Fig. [3] (in this section we consider a ACDM cosmology 
with (n m0 ,Q ao, a 8 ,n s ,h) = (0.27,0.73,0.79,0.95,0.7)). 
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Fig. 4. The Lagrangian map q' i— > r 1 given by the spheri- 
cal dynamics (neglecting shell-crossing) for the saddle-point 
(|2"5)) with a nonlinear density contrast 8 = 200 at the 
Eulerian radius r. We plot the curves obtained at z = 
for several masses, as in Fig. [3] Inner shells have already 
gone once through the center (hence their dynamics is no 
longer exactly given by Ea. lfTTj) ) but they have not crossed 
the radius r yet. 

We show in Fig. [4] the Lagrangian map, q' > r', given 
by the spherical dynamics (i.e. the function J-(Sl) where 
we neglect shell-crossing) for the saddle-point (|2"5")) . with a 
nonlinear density contrast S — 200 at the Eulerian radius 
r, at redshift = 0. Inner shells have already gone once 
through the center but they have not reached radius r yet. 
Even though their dynamics is no longer exactly given by 
Eq. (ITT1) . an exact computation would give the same prop- 
erty as the increasing mass seen by these particles, as they 
pass outer shells, should slow them down as compared with 
the constant-mass dynamics. In agreement with Fig. [31 for 
larger masses, which have a larger central linear density 
contrast, shell-crossing has moved to larger radii (the local 
maximum of r'/r, to the left of r' — 0, is higher). 

From the Lagrangian map, q' H ¥ r' , we define the nonlin- 
ear density threshold, S + , as the nonlinear density contrast 
S r reached within radius r at the time when inner shells 
first cross this radius r. Then, up to <5+, the mass within 
the Lagrangian shell q has remained constant, so that the 
saddle-point ([2"5j) - (|2"6")) is exact (this slightly underestimates 
<5+ as the expansion of inner shells should be somewhat 
slowed down by the mass of the outer shells they have over- 
taken). We show in Fig.[5]the dependence on redshift of this 
threshold S + , for several masses. In agreement with Figs. [31 
[H this threshold is smaller for larger masses. It shows a 
slight decrease at higher redshift as £l m (z) grows to unity. 
We can see that for massive clusters at z = 0, which have 
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M= 1 O 11 h - 1 M r 



M=10 l5 h-'M 



Fig. 5. The nonlinear density contrast 5+, beyond which 
shell-crossing must be taken into account, as a function of 
redshift for several mass scales. 

a mass of order 10 15 /i _1 Mq, the density contrast 6 ~ 200 
should separate outer shells with radial accretion from inner 
shells with a significant transverse velocity dispersion, built 
by the radial-orbit instability that dominates the dynamics 
after shell-crossing, see Valageas 2002b (in particular the 
appendix A). This agrees with numerical simulations, see 
for instance the lower panel of Fig. 3 in Cuesta et al. (2008), 
which show that rare massive clusters exhibit a strong ra- 
dial infall pattern, with a low velocity dispersion, beyond 
the virial radius (where S r ~ 200), while inner radii show 
a large velocity dispersion (even though we can distinguish 
close to the virial radius the outward velocity associated 
with the shells that have gone once through the center). 
Small-mass halos do not show such a clear infall pattern 
and the velocity dispersion is significant at all radii (e.g., 
upper panel of Fig. 3 in Cuesta et al. 2008). This expresses 
the fact that such halos, associated with typical events and 
moderate density fluctuations, are no longer governed by 
the spherical saddle-point ([2"5]l. Indeed, at low mass and 
small Lagrangian radius q, a q is no longer very small so 
that the path integral (|2"2")l is no longer dominated by its 
minimum and one must integrate over all typical initial 
conditions, including large deviations from spherical sym- 
metry. 

Next, the amplitude bh q and the minimum tp are given 
as functions of y by the implicit equations (Legendre trans- 
form) 



V 



-t^t and ip = yy + — 



dg 



(27) 



where the function r((?) is defined from the spherical dy- 
namics through the parametric system (Valageas 2002a; 
Bernardeau 1994b) 



g = S r = F{8Lq) and r = -5 Lq 
The system ((27)) also reads as 



ip = mm 

g 



yQ 



(28) 



(29) 



Note that the function t(G), whence the cumulant gener- 
ating function (p(y), depends on the shape of the linear 



power spectrum through the ratio of linear power a r jo q at 
Eulcrian and Lagrangian scales r and q. Finally, from the 
inverse Laplace transform (I2TJ1) . a second steepest-descent 
integration over the variable y gives (Valageas 2002a) 



cr r -> : V{S r ) ~ e 



-*L/(20 



(30) 



Thus, in the quasi-linear limit, the distribution V(S r ) is gov- 
erned at leading order by the Gaussian weight of the linear 
fluctuation 8^ q that is associated with the nonlinear density 
contrast S r through the spherical dynamics. We can note 
that Eq. (j3"0)) could also be obtained from a Lagrange multi- 
plier method, without introducing the generating function 
tp{y). Indeed, in the rare-event limit, where V{5) is gov- 
erned by a single (or a few) initial configuration, we may 
write 



rare events : V(S) 



max e 
{Mq]|Mfc,]=<5} 



\5l.C7\Sl 



(31) 



That is, V{5) is governed by the maximum of the Gaussian 
weight e -( <5 L' C £ 1 ' l5 i)/ 2 subject to the constraint 5 t [8l\ = 5 
(assuming there are no degenerate maxima). Then, we can 
obtain this maximum by minimizing the action S[Si]/a^ 
of Eg. (1231) . where y plays the role of a Lagrange multi- 
plier. This gives the saddle-point (|2"5j) . and the amplitude 
8L q and the radius q are directly obtained from the con- 
straint 5 = F{8Lq), as in Eq. (|2"o) . Then, we do not need 
the explicit expression of the Lagrange multiplier y, as this 
is sufficient to obtain the last expression of the asymptotic 
tail (|30l) . Nevertheless, it is useful to introduce the generat- 
ing function <p(y), which makes it clear that the Lagrange 
multiplier y is also the Laplace conjugate of the nonlinear 
density contrast S as in Eq. (fT9")) , since it is also of inter- 
est by itself, as it yields the density cumulants through the 
expansion (l2Tj) . 

The asymptotic (j30|) holds in the rare-event limit. This 
corresponds to both the quasi- linear limit, ay — > at fixed 
density contrast S r , and to the low-density limit, S r — > 
at fixed ay, as long as there is no shell-crossing. This lat- 
ter requirement gives a lower boundary <5_ , for linear power 
spectra with a slope n > — 1, and an upper boundary S + , for 
any linear spectrum (Valageas 2002b). This upper bound- 
ary was shown in Fig. [5] for a ACDM cosmology, for several 
masses. 

Indeed, at large positive density contrast, shell-crossing 
always occurs, as seen in Figs. |4][5] This invalidates the 
mapping Sl h> S = J~(8l) obtained from Eq. fTTl) . as mass 
is no longer conserved within the Lagrangian shell q, so 
that the asymptotic behavior ([30]) is no longer exact. In 
fact, as shown in Valageas (2002b), after shell-crossing it is 
no longer sufficient to follow the spherical dynamics, even if 
we take into account shell-crossing. Indeed, a strong radial- 
orbit instability develops so that the sensitivity to initial 
perturbations actually diverges when particles cross the 
center of the halo. Then, the functional d r [dj,(q)] is sin- 
gular at such spherical states (i.e. it is discontinuous as 
infinitesimal deviations from spherical symmetry lead to a 
finite change of 6 r ) and the path integral (|22|) is no longer 
governed by spherical states that have a zero measure. As 
noticed above, this also means that, in the limit of rare 
massive halos, the nonlinear density threshold 5+ separates 
outer shells with a smooth radial flow from inner shells with 
a significant transverse velocity dispersion. Thus, 5 + marks 
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the virialization radius where isotropization of the veloc- 
ity tensor becomes important, in agreement with numerical 
simulations (e.g., Cuesta et al. 2008). 

It is interesting to note that a similar approach can be 
developed for the "adhesion model" (Gurbatov et al. 1989), 
where particles move according to the Zcldovich dynamics 
(Zeldovich 1970) but do not cross because of an infinitesi- 
mal viscosity (i.e. they follow the Burgers dynamics in the 
inviscid limit). Moreover, in the one-dimensional case, with 
a linear power-spectrum slope n = —2 or n = (i.e. the 
linear velocity is a Brownian motion or a white noise), it 
is possible to derive the exact distribution V(S r ) by other 
techniques and to check that it agrees with the asymptotic 
tail (|5D|) . as seen in Valageas (2009a, b). 

4. Mass function of collapsed halos 

The quasi-linear limit (T50|) of the distribution V(S r ) clearly 
governs the large-mass tail of the mass function n(M)dM, 
where we define halos as spherical objects with a fixed den- 
sity contrast S r = S. Indeed, since the Lagrangian radius 
q is related to the halo mass by M = p m 47rg 3 /3, massive 
halos correspond to large Eulerian and Lagrangian radii, 
and the limit M — > oo corresponds to the quasi-linear limit 
ay — > 0. Then, going from the distribution V(S r ) to the 
mass function n(M) can introduce geometrical power-law 
prefactors since halos are not exactly centered on the cells 
of a fixed grid, as discussed in Betancort-Rijo & Montero- 
Dorta (2006), but the exponential cutoff remains the same 
as in Eg. ([50)1 . whence 



M oo : ln[n(M)] 



2a 2 



with 5 L =T- 1 (5), (32) 



where a — a q with M — p m 47r<7 3 /3. A simple approxi- 
mation for the mass function that satisfies this large-mass 
falloff can be obtained from the Press-Schechter method 
(PS), see Press & Schechter (1974), but using the actual 
linear threshold Sl = J r_1 (<5) associated with the nonlin- 
ear threshold S that defines the halo radius r, rather than 
the linear threshold S c ~ 1.686 associated with the spheri- 
cal collapse down to a point. This yields for the number of 
halos in the range [M, M + dM] per unit volume, 



n(M)dM=^f(a) 



da 



with 



/(*) 



2 S, 



7r a 



(33) 



(34) 



The mass function (|3~4"|) includes the usual prefactor 2 that 
gives the normalization 



da 



/(") = 1. 



(35) 



which ensures that all the mass is contained in such halos. 
However, we must note that the power-law prefactor in (|34[) 
has no strong theoretical justification since only the expo- 
nential cutoff ([3"2l was exactly derived from (|30l) . In princi- 
ple, it could be possible to derive subleading terms (whence 
power-law prefactors) in the quasi-linear limit for V(S r ), by 
expanding the path integral (|2"2")l around the saddle-point 
(|2"5"|) . Then, one would also need to take more care of the 



prefactors that would arise as one goes from the distribution 
V(S r ) to the mass function n(M). However, this expansion 
would not allow one to derive the shape of the mass func- 
tion at low masses, as this regime is far from the quasi-linear 
limit and involves multiple mergers far from spherical sym- 
metry. Until a new method is devised to handle this regime, 
one must treat the prefactors as free parameters to be fit- 
ted to numerical simulations, in order to describe the mass 
function from small to large masses. 

Here we can note that in the Press-Schechter approach 
(and in most models) the mass function ([34]) is obtained 
from the linear density reached within the sphere centered 
on each mass element. That is, a particle is assumed to be 
part of a halo of mass larger than M if the sphere of mass 
M (or a sphere of mass larger than M in the excursion set 
approach of Bond et al. 1991), centered on this particle, has 
a linear density contrast above a threshold Sl- As pointed 
out in Sheth, Mo & Tormen (2001), and discussed in their 
section 3, since all particles are not located at the cen- 
ter of their parent halo, the correct criteria should rather 
be that the particle belongs to a sphere, not necessarily 
centered on this point, of mass greater than M, that has 
collapsed by the redshift of interest. Within the excursion 
set approach of Bond et al. (1991), this means that one 
should consider the first-crossing distributions associated 
with center-of-mass particles, rather than with randomly 
chosen particles, as argued in Sheth, Mo & Tormen (2001). 
Then, the latter suggest that this could modify the numer- 
ical factor in the exponential tail of the mass function, that 
is, lead to a smaller factor Sl in Eq. (|3"2"|) than the one asso- 
ciated with spherical (or ellipsoidal) collapse dynamics. We 
point out here that this effect should only give subleading 
corrections to the tail (|3"2"|) . so that the factor Sl in Eg. (1521 
remains exactly given by the spherical collapse dynamics. 

This can be seen from the fact that the same effect 
would apply to the density probability distribution V(S r ), 
as randomly placed cells of radius r are typically not cen- 
tered on halo profiles. Nevertheless, the results (|2"T)) - (|3ip are 
exact in the quasi-linear limit (as can also be checked by 
the comparison with standard perturbation theory for the 
cumulant generating function <p(y)), and off-center effects 
would be included in the subleading terms, computed as 
usual by expanding the path integral (|22l) around its saddle- 
point, which would typically give the power-law prefactor to 
the tail (|3U| . In terms of the halo mass function itself, this 
point was also studied in Betancort-Rijo & Montero-Dorta 
(2006), who found as expected that at large mass such a 
geometrical factor only modifies that power-law prefactor. 

As for the probability distribution V(S r ), it is interesting 
to note that a similar high-mass tail can be derived for the 
"adhesion model" , where halos are defined as zero-size ob- 
jects (shocks). Again, in the one-dimensional case, for both 
ri = —2 and n = 0, where the exact mass function can be 
obtained by other means, one can check that it agrees with 
the analog of the asymptotic tail (|3"2")l (Valageas 2009a, b,c). 
Thus, we can check that in these two non-trivial examples 
the leading-order terms for the large-mass decays of V{S r ) 
and n(m) are exactly set by saddle-point properties and are 
not modified by the off-center effects discussed above. 

We can note that the explanation of the large-mass tail 
(f3"2"|) by the exact asymptotic result of the steepest-descent 
method described in section [3] agrees nicely with numerical 
simulations. This is most clearly seen in Figs. 3 and 6 of 
Robertson et al. (2009), who trace back the linear density 
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contrast Sl of the Lagrangian regions that form halos at 
z = 0. Their results show that the distribution of linear 
contrasts Sl, measured as a function of mass (or of a(M)), 
has a roughly constant lower bound 5^, with 5^ ~ 1.6, and 
an upper bound cT^that grows with a(M). We must note 
however that the difference between 1.59 and 1.6754 (which 
would be the standard threshold in their ACDM cosmology 
with Q m o — 0.27) is too small to discriminate both values 
from the results shown in their figures, so that these numer- 
ical simulations alone do not give the asymptotic value 5^ 
to better than about 0.2. They obtain the same results when 
they define halos by nonlinear density contrasts S — 100 or 
S = 600, with a lower bound S~[ that grows somewhat with 
5. In terms of the approach described above, this behavior 
expresses the fact that the most probable way to build a 
massive halo of nonlinear density contrast 6 is to start from 
a Lagrangian region of linear density contrast Sl = J-~ 1 (S), 
which obeys the spherical profile (|25p and corresponds to 
the saddle-point of the action (|2"5|) . As recalled in section [31 
the path integral (j2"2"|) is increasingly sharply peaked around 
this initial state at large mass scales, which explains why 
the dispersion of linear density contrasts Sl measured in the 
simulations decreases with a(M). At smaller mass, one is 
sensitive to an increasingly broad region around the saddle- 
point, which mostly includes non-spherical initial fluctua- 
tions. Since these initial conditions are less efficient to con- 
centrate matter in a small region one needs a larger linear 
density contrast Sl to reach the same nonlinear threshold 
S within the Eulerian radius r, which is why the distribu- 
tion is not symmetric and mostly broadens by increasing its 
typical upper bound 5^- In principles, it may be possible 
to estimate the width of this distribution, at large masses, 
by expanding the action S[Sl] around its saddle-point. 

We compare in Fig. [6] the prediction (J34J) (solid line 
labeled as Sl = J-^ 1 {S)) with results from numerical sim- 
ulations, for the nonlinear threshold 5 — 200 at rcdshift 
z = 0. In our case, this corresponds to a linear density 
contrast Sl — 1-59, that is obtained from Fig. [5] or the 
fit (fT2")) . As in section [21 we consider a ACDM cosmol- 
ogy with (n m(h n A o,<j 8 ,n s ,h) = (0.27,0.73,0.79,0.95,0.7). 
This corresponds to the cosmological parameters of the 
largest-box numerical simulations of Tinker et al. (2008), 
which allow the best comparison with the theoretical pre- 
dictions, as they also define halos by the density thresh- 
old S — 200 with a spherical-overdensity algorithm. The 
numerical results are the fits to the mass function given 
in Sheth & Tormen (1999) ("ST"), Jenkins et al. (2001) 
("J"), Reed et al. (2003) ( "R" ) , Warren et al. (2006) ("W") 
and Tinker et al. (2008) ("T"). Note that these mass func- 
tions are defined in slightly different fashions, using ei- 
ther a spherical-overdensity or friends-of-friends algorithm, 
and density contrast thresholds that vary somewhat about 
S = 200. However, they agree rather well, as the depen- 
dence on S is rather weak. This can be understood from 
Fig. [U which shows that around S — 200 the linear den- 
sity contrast Sl = J-^ 1 (S) has a very weak dependence 
on 5. We also plot the usual Press-Schechter prediction 
(dashed line labeled S c ), that amounts to replace Sl by 
S c = 1.6754 in Eq.® (since Cl m = 0.27 at z = 0). We 
can see that using the exact value Sl = J-^ 1 (S) signifi- 
cantly improves the agreement with numerical simulations 
at large masses. Note that there are no free parameters 
in Eg. ([3"4"| . Of course, at small masses the mass function 
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Fig. 6. The mass function at redshift z — of halos de- 
fined by the nonlinear density contrast S = 200. The points 
are the fits to numerical simulations from Sheth & Tormen 
(1999), Jenkins et al. (2001), Reed et al. (2003), Warren et 
al. (2006) and Tinker et al. (2008). The dashed line ("<5 C ") is 
the usual Press-Schechter mass function, while the solid line 
that goes close to the Press-Schechter prediction at small 
mass is the mass function (l3~4")l with the exact cutoff (|3"2")l . 
Here we have Sl — 1.59. The second solid line that agrees 
with simulations over the whole range is the fit (|3"6"1) . that 
obeys the same large-mass exponential cutoff. 

([34|) closely follows the usual Press-Schechter prediction 
and shows the same level of disagreement with numerical 
simulations. This is expected since only the exponential cut- 
off ([32]) has been exactly derived from section [3] Since at 
large masses the power-law prefactor 1 ja in Eq. (|3"4")l is also 
unlikely to be correct, as discussed above, we give a simple 
fit that matches the numerical simulations from small to 
large masses (solid line that runs through the simulation 
points) while keeping the exact exponential cutoff: 

f(a) = 0.5 [(0.6 is) 2 - 5 + (0.62 */) a5 ] e^/ 2 , (36) 
with 

v = — , S L = T~ x {S) ~ 1.59 for S = 200, Q m = 0.27(37) 
a 

At higher redshift or for other f2 m one simply needs to use 
the relevant mapping T~ (S), shown in Fig. [2]or given by 
the simple fit (|12l) for S = 200. This mass function also 
satisfies the normalization (|35[) . whatever the value of the 
threshold Sl- 

Thus, we suggest that mass functions of virialized halos 
should be defined with a fixed nonlinear density threshold, 
such as S = 200, and fits to numerical simulations should 
use the exact exponential cutoff of Eq. (|3"2"]) . with the ap- 
propriate linear density contrast Sl = J'^ 1 (S), rather than 
treating this as a free parameter. This would automatically 
ensure that the large mass tail has the right form (up to sub- 
leading terms such as power-law prefactors), as emphasized 
by the reasonable agreement with simulations of Eq. (j34|) at 
large masses. Moreover, it is best no to introduce unneces- 
sary free parameters that become partly degenerate. Here 
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Fig. 7. The mass functions at z = of halos defined by 
the nonlinear density contrasts S — 100 (upper panel) and 
5 = 300 (lower panel) . The points show the numerical sim- 
ulations of Tinker et al. (2008). The dashed line is the usual 
Press-Schechter mass function, while the solid lines corre- 
spond to Eqs.® and fl35J|, with now S L = J 7 ' 1 (100) ~ 
1.55 (upper panel) and Sl = J-~ 1 (300) ~ 1.61 (lower 
panel). 



we must note that Barkana (2004) had already noticed that, 
taking the spherical collapse at face value and defining ha- 
los by a nonlinear threshold 6 (he chose S — 187r 2 ), one 
should use the linear threshold Sl = J-~ 1 (6) for the Press- 
Schechter mass function. Unfortunately, noticing that the 
value of 6l given by fits to numerical simulations was even 
lower, he concluded that this was not sufficient to reconcile 
theoretical predictions with numerical results. As shown by 
Fig. O this is not the case, as the parameter Sl used in 
fitting formulae is partly degenerate with the exponents of 
the power-law prefactors, so that it is possible to match 
numerical simulations while satisfying the large-mass tail 
l|32p. Of course, the actual justification of the asymptote 
l|32p is provided by the analysis of section [3[ which shows 
that below the upper bound <5+ the spherical collapse is in- 
deed relevant and asymptotically correct at large masses, 
as it corresponds to the saddle-point of the action (f2"3"f . 

Note that the result (|32p also implies that the mass 
function f(a) is not exactly universal, since the mapping 
S i y Sl ~~ J-^ 1 (S) shows a (very) weak dependence on 
cosmological parameters, see Fig. [5] Using the exact tail 



([3"2"jl should also improve the robustness of the mass func- 
tion with respect to changes of cosmological parameters and 
redshifts. 

Next, we compare in Fig.[7]the mass functions obtained 
at z — for the density thresholds S = 100 and S — 300 
with the numerical simulations from Tinker et al. (2008), 
who also considered these density thresholds. We show the 
usual Press-Schechter mass function (dashed line) and the 
results obtained from Eqs. (|3"4")) and ([551) (solid lines), with 
now S L = ^(lOO) ~ 1.55 and S L = J 7-1 (300) ~ 1.61. 
We can see that the large mass tail remains consistent with 
the simulations, but for the case 5 = 100 it seems that the 
shape of the mass function at low and intermediate masses 
is modified and cannot be absorbed through the rescaling of 
8l- It appears to follow Eq. (|3"4"|) rather than the fit (|3l)]) . but 
this is likely to be a mere coincidence. We should note that 
Fig. [5] shows that 5+ < 300 for massive halos, so that shell- 
crossing should be taken into account and the tail (f3"2"j) is 
no longer exact for S — 300, although it should still provide 
a reasonable approximation, as checked in Fig. [JJ Thus, 
even though Tinker et al. (2008) also studied higher density 
thresholds, we do not consider such cases here as the tail 
(1321) no longer applies. 



5. Halo density profile 

As for the mass function n(M), the analysis of section [3] 
shows that the density profile of rare massive halos is given 
by the spherical saddle-point (f2"5j) . see also Barkana (2004) 
and Prada et al. (2006). This holds for halos selected by 
some nonlinear density threshold S in the limit of rare 
events, provided shell-crossing has not occurred beyond the 
associated radius r. In particular, this only applies to the 
outer part of the halo since in the inner part, at r' <C r, 
shell-crossing must be taken into account. Then, as dis- 
cussed in section [3[ a strong radial-orbit instability comes 
into play and modifies the profile in this inner region, as 
deviations from spherical symmetry govern the dynamics 
and the virialization process (Valageas 2002b). 

We compare in Fig. \8\ the nonlinear density profile ob- 
tained from Eq. (l25[) with fits to numerical simulations. 
We plot the overdensity within radius r' , 1 + S r = p m (< 
r')/p m , as a function of radius r. This is again obtained 
from Eg. ([2"5]l with the mapping q 13 = (1 + S r r)r' 3 and 
5 r i = J-(SLq')- Since this only applies to the limit of rare 
events, we choose for each mass M a redshift z such that 
cr(M, z) = 0.5. Thus, smaller masses are associated with 
higher redshifts. The results do not significantly depend on 
the precise value of the criterium used to define rare events, 
here cr(M, z) = 0.5. The points in Fig.[5]are the results ob- 
tained from a Navarro, Frenk & White profile (Navarro et 
al. 1997, NFW), 



p(r') 



(r'/r^l + r'/r^' 
or an Einasto profile (Einasto 1965), 



Hp(r')/ P - 2 } = — [(r7r_ a r-l] 
a 



(38) 



(39) 



For the NFW profile, the characteristic radius r s is ob- 
tained from the concentration parameter, c(M,z) — r/r s , 
where as in the previous sections r is the Eulerian ra- 
dius where S r — 200 (i.e. r = r2oo)- Then, we use in 
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Fig. 8. The density profile of rare massive halos, for sev- 
eral masses M. For each mass the redshift is such that 
(t(M, z) — 0.5, as the theoretical prediction from Eq. (j23)) 
(solid line) only applies to rare events. The second branch 
that appears at small radii is due to the shell-crossing, hence 
the theoretical prediction only holds to the right of this 
branch. The points are fits to numerical simulations, based 
on an NFW profile (Dolag et al. 2004; Duffy et al. 2008) or 
an Einasto profile (Duffy et al. 2008). Note that we show 
the mean density within radius r', p m (< r ')/Pnu rather 
than the local density at radius r' , so that the NFW and 
Einasto profiles are integrated once. 



Fig. [5] the two fits obtained in Dolag et al. (2004) and 
Duffy et al. (2008) from numerical simulations, of the form 
c(M, z) = a(M/M a ) b (l + z) c . For the Einasto profile, we use 
the fits obtained in Duffy et al. (2008) for the concentration 
parameter, c(M, z) = r/r_2, and for the exponent a (writ- 
ten as a quadratic polynomial over v = S c /a(M, z) as in 
Gao et al. (2008)). Then, we integrate over r' the profiles 
to obtain the overdensity profiles p m (< r')/p m 
shown in Fig. [8] (as p m (< f') is the mean density within ra- 
dius r' and not the local density at radius r' as in Eqs. ([35|) - 

62)). 

We can check in Fig. [5] that our results agree reasonably 
well with these fits to numerical simulations over the range 
where both predictions are valid. Note that our prediction 
has no free parameter, since it is given by the saddle-point 
profile (1231) . At large radii we recover the mean density of 
the Universe, while the numerical profiles (|38p - (j39|) go to 
zero, but this is an artefact of the forms (|38l) - (j3"9")l . that 
were designed to give a sharp boundary for the halos and 
are mostly used for the high-density regions. For a study of 
numerical simulations at outer radii, see Figs. 8 and 10 of 
Prada et al. (2006) which recover the mean density of the 
Universe at large scales. At small radii, the second branch 
that makes a turn somewhat below r is due to shell-crossing 
that makes the function p m (< r 1 ) bivaluate. Below the max- 



imum shell-crossing radius, for each Eulerian radius r' there 
are two Lagrangian radii q' , a large one that corresponds 
to shells that are still falling in, and a smaller one that cor- 
responds to shells that have already gone once through the 
center (note that in agreement with Fig. [J] shell-crossing 
appears at slightly smaller relative radii for smaller mass). 
Then, the theoretical prediction only holds to the right of 
this second branch, where there is only one branch and 
no shell-crossing. Note that the theoretical prediction (1231) 
explicitly shows that the halo density profiles are not uni- 
versal. Within the phenomenological fits (j3"51) - (j3"9")) this is 
parameterized through the dependence on mass and red- 
shift of the concentration parameter and of the exponent 
a. However, we can see that over the regime where Eg. (1231) 
applies the local slope of the halo density is p(r') ~ r'~ 2 , 
which explains the validity of the fits (T3"8l - (j3"9")) in this do- 
main. Unfortunately, our approach cannot shed light on the 
inner density profile, where 5 r i 3> 200, which is the region of 
interest for most practical applications of the fits (|38D - (j3"9")l . 

The same approach, based on the spherical collapse, was 
studied in greater details in Betancort-Rijo et al. (2006) 
and Prada et al. (2006). They consider the "typical" pro- 
file ([23| that we study here, as well as "mean" and "most 
probable" profiles. In agreement with the steepest-descent 
approach of section^ they find that the most probable pro- 
file closely follows the typical profile (123)) and they obtain a 
good match with numerical simulations for massive halos, 
paying particular attention to outer radii. Therefore, we do 
not further discuss halo profiles here. 



6. Halo bias 

In addition to their multiplicity and their density profile, a 
key property of virialized halos is their two-point correla- 
tion function. At large scales it is usually proportional to 
the matter density correlation, up to a multiplicative fac- 
tor b 2 , called the bias of the specific halo population. We 
revisit in this section the derivation of the bias of massive 
halos, following Kaiser (1984), and we point out that pay- 
ing attention to some details it is possible to reconcile the 
theoretical predictions with numerical simulations, without 
introducing any free parameter. 

As seen in the previous sections, rare massive halos can 
be identified with rare spherical fluctuations in the initial 
(linear) density field. More precisely, as in section[3]we may 
consider the bivariate density distribution, V(5 ri , <5 r2 ), of 
the density contrasts 5 ri ,S r2 , in the spheres of radii ri,r 2 , 
centered at points Xi,x 2 . Thus, we introduce as in Ea. (IT9l) 
the double Laplace transform </?(yi, 2/2), 

e -v(yi,y2)/<? 2 _ ^ p -(yi$r 1 +y25r 2 )/°- 2 ^ 

dS ri dS r2 e ^ lS ^ +v^)/^ 2 p(S ri ,5 r2 ), (40) 



where a 1 = a 2 r (0) is the linear variance flT8]) at some scale 
r. If r*i = r-i we may take r = T\ — r2, otherwise it can 
be any intermediate scale, as the results do not depend on 
this factor. The only requirement is that a should scale 
as cr ri and o> 2 , which is obtained by choosing for instance 
a fixed ratio r/^/rirJ. Then, the probability distribution 
V(S ri , S r2 ) can be obtained from the cumulant generating 
function ^(2/1,2/2) through a double inverse Laplace trans- 
form, as in Ea. (|20p . Next, <p(j/i, 2/2) can be written in terms 
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of the linear density field as a path integral, such as (12"2"1) . 
with an action that now reads as 



S[S L ] = yi 5 ri [S L ] + y 2 S r2 [S L ] + y 6 l-C l \S l 



(41) 



Then, for rare events the tail of the distribution V(8 ri , 8 r2 ) 
is governed by the last Gaussian weight of (|4T|) . as in 
Eq.®, 



0: V(S ri ,S r2 ) 



-h$L C , 1 .S L 



(42) 



where <5i[q] is the relevant saddle-point of the action S. 
Unfortunately, since the action (|4"Tj) is no longer spheri- 
cally symmetric, it is not possible to obtain an explicit 
expression of its minimum. Therefore, we must rely on 
some simple approximations. In the limit of large distance, 
%12 = |x2 — xx| — > oo, between the positions of both 
halos, the linear field #l(s) around the Lagrangian po- 
sitions Si,S2, of the halos should follow the profile (|2"5|) 
(we note the Lagrangian coordinates to avoid confusion 
with the Lagrangian radii qi). Thus, we neglect the tidal 
forces of the halos, but we keep track of the mean displace- 
ment of each halo, due to the gravitational attraction from 
the other one, by distinguishing their Lagrangian positions 
Si, S2, from their Eulerian positions Xi, X2. Within this ap- 
proximation, the distribution V(8 ri , 8 T2 ) can be estimated 
from the Gaussian distribution Vh(8h qi ^8h q2 ) of the lin- 
ear density contrasts 8L qi ,SL q2 , at positions Si,S2, within 
the Lagrangian spheres of radii q\ and (72 , with the mapping 
(j2"6")l . This closely follows the approach introduced in Kaiser 
(1984), except for the distinction between and x^. The 
joint distribution Vhi^Lq^^hq-i) reads as (see also Kaiser 
1984; Politzer & Wise 1984) 



VL{5 Lqil 8 L q 



1 



Ar 



(2ir)VdetM 
where M is the linear covariance matrix, 



M 



'12 



2 

^12 



(43) 



(44) 



Here we defined from Eq. (fT8"l) . of = a q . qi (0) and a\ 2 = 
a~ 1 q2 (s\ 2 ) with s%2 = |s2 — Si|. The inverse of the matrix 
M writes as 



'12 



2 „2 



err a 



'12 



_2 
"°12 



(45) 



This yields 

V L {S Lqi ,5 Lq2 ) = V L {SLq 1 )'PL{SLq 2 ) 

2S Lqi SLq 2 crh~Sl qi at 2 /af 



CT1O2 



'12 



x exp 



4) 



(46) 



Next, defining the real-space halo correlation ^u 1 ,M 2 ( r ) as 
the fractional excess of halo pairs (Kaiser 1984; Peebles 
1980), 

«x 1 ,x 2 (Mi,M 2 )dMidA/ 2 dxidx2 = 

n{M 1 )n(M 2 ) (1 + £ Ml ,M 2 (^12)) dAfidM 2 dxidx 2 , (47) 



which also gives the conditional probability, 
n Xl ,x 2 (M 2 |Mi) = n(M 2 ) (1 + ^(m)) , 



we write 



1 +€m 1 ,m 3 (xi2) ~ + 8 M (x 12 )) 



Vl(8 l1 ,6 L2 ) 
V L {8 L1 )V L {8 L2 y 



(48) 



(49) 



where the mass of each halo is given by Mj = p m 47rg|/3 
and 5 Li = -7-" 1 (<5j ) are the linear density contrasts which 
are associated with the nonlinear density contrasts Si that 
define the halos, as in section [4] and Fig. In Ea. (l4T)l) 
the factor (1 + 8m(xi2)) models the effects associated 
with the mapping from Lagrangian to Eulerian space. 
Indeed, in the limit of rare massive halos and large separa- 
tion, the number of neighbors is conserved and we write 
TiEuia^ckc^ = "LagSi2dsi2 (for an integral form of the 
conservation of pairs, or neighbors, see Peebles 1980). We 
take dsi2 = (1 + 8m )dxi2, where we define 8m{x\2) as the 
local nonlinear density contrast at Eulerian distance x\ 2 
from a halo of mass M, obtained from the profile (1231) . and 
we choose M — max(Mi,M 2 ) (which obeys the symme- 
try Mi -H- M 2 ), whence q = m&x(qi,q 2 ). This reflects the 
fact that the gravitational attraction of a massive halo pulls 
matter towards it center with a strength that depends on 
distance, so that the Jacobian |dq/<9x| = (1 + 8) is differ- 
ent from unity. In a sense this is a tidal effect, as a locally 
volume-preserving displacement would not affect the num- 
ber densities, which is why we take for 8m (#12) the den- 
sity contrast at distance xi 2 and not the density contrast 
within radius 212. We may note that a local bias model 
(Mo & White 1996), using a peak-background split argu- 
ment (Efstathiou et al. 1988), would rather give a quadratic 
factor (l + 8i)(l + 8 2 ). However, we prefer to keep the linear 
factor (|4TJ)l as it also remains consistent in case of large dis- 
placementfl This also expresses the fact that all pairs can- 
not simultaneously get closer (as fluctuations grow some ob- 
jects move closer but they also further separate from other 
emerging groups). Then, defining the halo bias as 



^Mi,M 2 ( r ) 



£mi,m 2 (Q 



(50) 



where r = xi 2 and £(r) is the nonlinear matter correlation, 
we obtain the bias from Eq. (|4"9")l . Note that the bias (|5U| 
does not factorize, b M M 2 ( r ) ^ ^Afi ( r ")^M 2 (>")> because of 
the terms 012 and 8m- 

At large separation, r — > 00, we are in the linear regime, 
so that the matter correlation reads as £(r) ~ tToo( r )> an d 
the local density contrast 8m (t) is small and close to the 
linear density contrast 8l(s). Note that this is the density 
contrast at Lagrangian radius s, and not the mean density 
contrast within radius s. Therefore, it is related to Ea. ([25|) 
by 



^ 2 5 L {q) = #- 



,q 3 S Lq ) 



(51) 



2 For instance, let us consider a large system which can be 
subdivided into cells of two classes, {+, — }, with matter densi- 
ties {p+,p_} and volume fractions {r/+,ri-}, and p+ > p > P-, 
77+ < rj-. Then, from the conservation of volume (rj+ +rj- = 1) 
and mass (?7+p+ + Tj—p- = p), we obtain in the limits p+ ^> 1 
and p- « 1, £ = (p 2 )c/p 2 ~ p+/p- 
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whence 



Sl(s) = S L 



d 2 (q, s) 



(52) 



since at large distance we have 5 r > ~ &Lq' <C 1. Therefore, 
we obtain the Lagrangian separation s as the solution of 
the implicit equation 



with Sl = J- 1 (200) for instance, and using Ea.([T8"|) 



d- 2 {q,s) =4tt / dkk 2 P L (k)W(kq) 



sin(fcs) 



al (s). (53) 



Then, from Eq.([5 
same linear threshold 



& M 1 ,M 2 ( r ) 



1 



the bias of halos defined by the 
, = J r_1 (<5), reads as 

« -<, 92 ( S ) 



For equal-mass halos this simplifies as 



92 



(54) 



i 




x exp 



L u q,q 



2 W-<51< 9 M/^ 2 



9,9 



(s) 



(55) 



Note that the argument of the exponential is not necessar- 
ily small, as stressed in Politzer & Wise (1984). Indeed, we 
only assumed a large separation limit, i.e. a 2 , q {s) <C <r 2 , 
and a rare-event limit, Sh/vq 1- Thus, at fixed (small) 
ratio <j 2 q {s)/<j 2 , we obtain a large exponent in the limit 
of large-mass halos, a 2 — > 0. Then, keeping the full ex- 
pressions ([54]) or (1551) gives a nonlinear bias, since b 2 M (r) 
is not a simple number and shows a non-trivial scale de- 
pendence, as will be clearly seen in Fig. [TT] below. Indeed, 
these results cannot be recovered through a linear biasing 
scheme, where the fluctuations of the halo number density 
field, <5haio = (n — n)/n, are written as 5h a io{M) = b(M)Sji 
(where the matter density field is smoothed over some 
larger scale R), which would lead to £,Ai(r) oc o 2 R R (r) 
whereas Eos.([54 ]) - ([55"]) generate powers of all orders over 
<j 2 q (r). In the local bias framework (Fry & Gaztanaga 1993; 

Mo et al. 1997), where one writes <5 ha io(M) = Ei h(M)S R , 
this could be interpreted as non-zero bias coefficients bi for 
i > 2. However, Eqs. ([5"4"]) - ([5"5"]) are not equivalent to such a 
local model, inspired from a peak-background split argu- 
ment (Efstathiou et al. 1988). Indeed, they do not involve 
any external smoothing scale R and they depend on the 
three linear correlations <7q (r) , a 2 (r) and <J 2 q (r). 

Finally, we must express the Lagrangian separation s 
in terms of the Eulerian distance r. At lowest order we 
again consider each halo as a test particle that falls into 
the potential well built by the other halo (i.e. we neglect 
backreaction effects) . Then, from the analysis of section [3] 
and the linear profile ((25]) . we know that a test particle at 
Lagrangian distance q' from one halo has moved to position 
r', according to the mapping (l2l)|) . Using Eg. (1231) this gives 
at first order 



i 



5hq> 

3 



(56) 



r = s 1 — 



qua 



q2,s 



'12 



(57) 



where again we only kept the first-order term and we took 
into account the displacements of both halos. Together with 
Eq . (f5~4"]) . or Eq. ([55|) . this defines our prediction for the bias 
of massive collapsed halos. Note that this approach also 
applies to the cross-correlation between different rcdshifts. 

At large separation, r — > oo, and fixed mass (i.e. fixed 
Oq), we may linearize the bias (|55|) over o 2 {s) as 



b 2 M (r) 



'0,0 



(r) 



<?q,0(s) 



2 

c2 a q,q 



(58) 



For large masses, where a q -C 1, but not too large, so that 
the exponent in (|55p is still small, this gives 



r — > oo, M — > oo 



Mr) - % aqM 



<? 2 q co,o(0 



(59) 



Thus we recover the result of Kaiser (1984) and Mo 
& White (1996), except for the multiplicative factor 
a q,q{ s )/ (T o,o{ r )- It expresses the facts that halos only probe 
the linear density field smoothed over the Lagrangian scale 
q (i.e. the formation of a halo does not depend on wave- 
lengths much smaller than its radius) and that halos have 
moved from distance s to r by their mutual gravitational 
attraction. This factor yields a weak scale dependence for 
b(M), as the slope of the linear power spectrum slowly 
varies with scale r. We can also note that at lowest order 
over <7 2 (s) we may write the solution of Eq. (|57[) as 




92, r 



(60) 



which provides a simple explicit expression for s. 

We compare in Figs. [MTTI the bias obtained from 
Eqs. (15"5j) . (|57|) . with fits to numerical simulations. We first 
show in Fig. [9] our prediction as a function of halo mass, at 
redshift z — and distance r = 50/i _1 Mpc. This is typi- 
cally the scale that is considered in numerical simulations 
to compute the large-scale bias, as b(r) is expected to be 
almost constant at large scales (Kaiser 1984; Mo & White 
1996), see also Ea.([5^|. We also plot the standard theoret- 
ical prediction from Mo & White (1996) (dashed line). We 
can see that our result (|551) - (|57I) agrees rather well with 
numerical simulations and the popular fit from Sheth, Mo 
& Tormen (2001). As expected, it follows the trend of the 
prediction from Mo & White (1996), since both derivations 
follow the spirit of Kaiser (1984) (i.e. one identifies halos 
from overdensities in the linear density field) and they agree 
at large scale for rare massive halos, up to a factor of order 
unity, as seen in Eq. (|5"9")) . Note that Eqs. ([55 ]) - ([57 ]) only ap- 
ply to the rare-event limit, as for small objects the approx- 
imations used in the derivation no longer apply. In particu- 
lar, halos can no longer be considered as spherical isolated 
objects, and one should take into account merging effects. 
Note that this caveat also applies to other analytical ap- 
proaches, such as Kaiser (1984) and Mo & White (1996). 
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Fig. 9. The halo bias b(a) , as a function of a(M) , at redshift 
z = and distance r = 50/i -1 Mpc. The solid line is the 
theoretical prediction (|55|) - ([57|) . and the dashed line is the 
bias obtained in Mo & White (1996). The points are the 
fits to numerical simulations, from Sheth, Mo & Tormen 
(2001) (crosses) and Pillepich (2009) (circles). 
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Fig. 10. The halo bias b(r), as a function of scale r, at 
redshift z = and for several masses. The solid line is the 
theoretical prediction (|55|) - ([57|) . The crosses show the large- 
scale fit to numerical simulations from Sheth, Mo & Tormen 
(2001), while the triangles show the fit from Hamana et al. 
(2001). The dashed line is the linearized bias (f55|). the dot- 
dashed line is the nonlinear bias (|55|) where we set s = r 
while the dotted line uses Eq.® (for M = 10 15 /i _1 M Q in 
the three cases). We only plot our predictions for r > 2q. 



Then, our prediction should only be used for large masses, 
for instance such that b > 1. 

Next, we compare in Fig. [TU] the dependence on r of the 
bias ([55l) - ((57|) with the fit to numerical simulations from 
Hamana et al. (2001) (using their cosmological parameters) . 
The crosses are the large-scale limit given by Sheth, Mo & 



Fig. 11. The halo bias b(r), as a function of scale r, at 
redshift z = 10 and for several masses. As in Fig. [TUJ the 
solid line is the prediction (l55|) - (l57l) . while the crosses show 
the large-scale fit from Sheth, Mo & Tormen (2001), but 
the squares now show the fit to numerical simulations from 
Reed et al. (2009). The dashed line is the linearized bias 
([55]) . the dot-dashed line is the nonlinear bias (|55j) where 
we set s — r while the dotted line uses Ea. ljoTfl) (for M = 
10 1:l /i _1 Mq in the three cases). Again we only plot our 
predictions for r > 2q. 



Tormen (2001) and we only plot our prediction (solid lines) 
down to scale r — 2q, since it should only apply to large halo 
separations. We can see that the scale-dependence that we 
obtain is opposite to the one observed in the simulations. 
However, both are very weak and the prediction (|55|) - ([57| 
may still lie within error bars of numerical results. For the 
largest mass, M = 10 15 /i _1 M©, we also plot for illustration 
the linearized bias (|58p (lower dashed line), and the non- 
linear bias (|55p where we set s — r (upper dot-dashed line) 
or we use Eq. (|6T))) (lower dotted line). As expected, at very 
large scales the linearized bias ([55)1 agrees with the nonlin- 
ear expression (|55|) . Using the simpler Eq. ((6"0"]) also gives the 
same results at large scales, but the Lagrangian to Eulerian 
mapping still gives a non-negligible correction as shown by 
the upper dot-dashed line where we set s = r. In any case, 
it is always best to use the full expression (|55p - ([57|) . 

We compare in Fig. [11] our results at high redshift, 
z = 10, with the fit to numerical simulations from Reed 
et al. (2009) (using their cosmological parameters). Again, 
the crosses are the large-scale limit of Sheth, Mo & Tormen 
(2001) and we only plot our prediction (solid lines) down 
to scale r = 2q. For M = 10 11 /i~ 1 M© we also plot the 
linearized bias (|58[) (lower dashed line), and the nonlinear 
bias (|55l) where we set s = r (upper dot-dashed line) or we 
use Eq.® (dotted line). As noticed in Reed et al. (2009), 
the scale-dependence is much steeper than the one found 
at small redshifts and it is not consistent with the fits ob- 
tained at low z in Hamana et al. (2001) or Diaferio et al. 
(2003). This was interpreted as a breakdown of universal- 
ity for massive halos at high redshift by Reed et al. (2009). 
However, we can see that our prediction (|55)) - ([57| agrees 
reasonably well with their numerical results. Therefore, the 
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change of behavior of the bias b(r) between the two regimes 
studied in Figs. [10] and [TTJ can be understood from the stan- 
dard picture of massive halos arising from rare overdensities 
in the initial (linear) Gaussian density field, by using the 
same theoretical prediction (|55|) - ([57| that applies to any z. 
We can note that the linearized bias ([55]). or the approx- 
imation s = r, show strong deviations in this regime and 
disagree with the simulations. Therefore, one should use 
the nonlinear bias (|55|) - ([57|) (but using Eg. (1501) gives sim- 
ilar results) and one cannot neglect the correction due to 
the Lagrangian to Eulerian mapping that is associated with 
s H> r. 
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Fig. 12. The ratio b 2 (M u M 2 )/[b(M 1 )b(M 2 )], from 
Eqs. (l54"]) and (|55|) . at fixed geometrical mean 
M = \JM\M 2 . Solid lines are for distance and red- 
shift (r, z) = (50/i _1 Mpc, 0), whereas dashed lines are for 
(r, z) — (3/i -1 Mpc, 10). Curves are labelled by their fixed 
geometrical mean M. 



On the other hand, the comparison with the result from 
(|58|) shows that the steep dependence on scale is due to the 
nonlinear term in (|55[) . i.e. keeping the exponential factor. 
Indeed, as noticed below Ea. ([55|) . and in Politzer & Wise 
(1984), the derivation of the bias presented above only as- 
sumes a large separation between very massive (rare) halos, 
that is, cr 2 q (s) <C <r 2 and 8hjo q 3> 1. Therefore, for suffi- 
ciently massive objects (that correspond to a large bias), 
the variance a 2 can be small enough to make the expo- 
nent in Eq. (j5"5"]) of order unity or larger. Then, one needs to 
keep the nonlinear expression (|55[) rather than expanding 
the exponential as in Eg. ([55]). As noticed below Eq.(l55l). 
this implies a nonlinear biasing scheme as the halo cor- 
relation is not proportional to the matter correlation but 
shows a steeper scale dependence. Within the local bias 
framework (Fry & Gaztanaga 1993; Mo et al. 1997), this 
could be interpreted as non-zero higher order bias parame- 
ters bi in the expansion <5haio(-W) = Yli &i(-^0^k- However, 
the bias (|5"5)) cannot be exactly reduced to such a model 
(but one could certainly derive within such a framework a 
good approximation to Eq. ([55|) . restricted to some larger 
scale R, by using Ea. lpO"]) . writing the correlations a q0 (r) 



and a 2 q (r) in terms of tTQ (r), expanding over (TQ (r) and 
finally smoothing over the external scale R of interest). 

We should stress here that our prediction (|55p - ([57|) has 
no free parameter. This can lead to a slightly larger inac- 
curacy as compared with fits to simulations in the regime 
where the latter have been tested, as in Fig. [TUJ but this 
improves the robustness of the predictions as one consider 
other regimes (e.g. other cosmological parameters or other 
redshifts as in Fig. ITT]) . Therefore, we think that Eas. ([55l) - 
(|57p could provide a useful alternative to current fitting 
formulae, as they can be readily applied to any set of cos- 
mological parameters or redshifts. 

Finally, we show in Fig. [T5] the bias ratio 
b 2 (M 1 ,M 2 )/[b(M 1 )b(M 2 )}, as a function of the mass 
ratio M 2 /Mi, at fixed geometrical mean M = \J M\ M 2 . 
We consider several mass scales M, at distance 
and redshift (r, z) — (50/i _1 Mpc, 0) (solid lines) 
and (r,z) = (3/i -1 Mpc, 10) (dashed lines). This 
shows that making the factorized approximation, 
b 2 (Mi,M 2 ) ~ b(Mi)b(M 2 ), can lead to an error of 
up to 50% for a mass ratio M 2 jM\ ~ 100. Therefore, it is 
best to use the full Ea.([53J). 

7. Conclusion 

We have pointed out in this article that the large-mass 
exponential tail of the mass function of collapsed halos 
is exactly known, provided halos are defined as spherical 
overdensities above a nonlinear density contrast threshold 
8 (i.e. using a spherical overdensity algorithm in terms of 
numerical simulations). This arises from the fact that mas- 
sive rare events are governed by (almost) spherical fluctu- 
ations in the initial (linear) Gaussian density field (if one 
does not explicitly breaks statistical isotropy by looking for 
non-spherical quantities). This is most easily seen from a 
steepest-descent approach, which becomes asymptotically 
exact in the large-scale limit, applied to the action 5 as- 
sociated with the probability distribution of the nonlinear 
density contrast within spherical cells. This result holds for 
any nonlinear threshold S used to define halos, provided it 
is below an upper bound <5+ that marks the point where 
shell-crossing comes into play. For a standard ACDM cos- 
mology, 5+ typically grows from 200 to 600 as one goes from 
lO 15 ft _1 M0 to 10 11 hr 1 Mq (which also corresponds to in- 
creasing redshift). This dependence on mass is due to the 
change of slope of the linear power spectrum with scale. 

We have also noted that in two similar systems, the one- 
dimensional adhesion models with Brownian or white-noise 
initial (linear) velocity, the same method can be used for 
both the density distribution V(8 r ) and the mass function 
n(M), and one can check that this yields rare-event tails 
that agree with the exact distributions, which can be de- 
rived by other techniques (Valageas 2009a, b). 

Therefore, defining collapsed halos by a threshold 8 = 
200 to follow the common practice, the large-mass tail 
of the halo mass function is of the form e^hl\ 2a ( M )) ) 
up to subleading prefactors such as power laws, where 
6l = J-^ 1 (S) is the linear density contrast associated to 8 
through the spherical collapse dynamics. In particular, we 
obtain 8l — 1-59 for 8 — 200. We checked that this value, 
which is slightly lower than the commonly used value of 
8 C = 1.686 associated with complete collapse, gives a good 
match with numerical simulations (at large masses) when 
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we simply use the Press-Schechter functional form. We also 
give a fitting formula, which obeys this exact exponential 
cutoff, that agrees with simulations over all mass scales. 

We suggest that halos should be defined by such a non- 
linear density threshold (i.e. friends-of-friends algorithms 
are not so clearly related to theoretical computations) and 
fits to numerical simulations should use this exact exponen- 
tial tail, rather than treating Sl as a free parameter. This 
would avoid introducing unnecessary degeneracies between 
fitting parameters and it would make the fits more robust. 

Next, we have briefly recalled that in the large- mass 
limit the outer density profile of collapsed halos is given 
by the radial profile of the relevant spherical saddle-point. 
This applies to radii beyond the density threshold 6+ , where 
shell-crossing comes into play. In agreement with numeri- 
cal simulations, for rare massive halos this separates an 
outer region dominated by a radial flow from an inner re- 
gion where virialization takes place and a strong transverse 
velocity dispersion quickly builds up. We have recalled that 
this can be explained from a strong radial-orbit instability, 
which implies that infinitesimal deviations from spherical 
symmetry are sufficient to govern the dynamics. 

Finally, following the approach of Kaiser (1984), we have 
obtained an analytical formula for the bias of massive ha- 
los that improves the match with numerical simulations. 
In particular, it captures the steepening of the scale de- 
pendence that is observed for large-mass halos at higher 
rcdshifts. This requires keeping the bias in its nonlinear 
form and taking care of the Lagrangian-Eulerian mapping. 
We also note that using a factorization approximation, 
b 2 (M 1 ,M 2 ) ~ b(M 1 )b(M 2 ), may lead to non-negligible in- 
accuracies. We stress that this analytical estimate of the 
bias contains no free parameter. Although this can yield 
a match to numerical simulations that is not as good as 
fitting formulae derived from the same set of simulations, 
it provides a more robust prediction for general cases, as 
shown by the good agreement obtained at both low and 
high redshifts (z = and z — 10), whereas published fitting 
formulae cannot reproduce both cases. We think this makes 
such a model useful for cosmological purposes, where it is 
desirable to have versatile analytical estimates that follow 
the correct trends as one varies cosmological parameters 
or redshifts. In particular, the scale-dependence of the bias 
of massive halos has recently been proposed as a test of 
primordial non-Gaussianity (e.g., Dalai et al. 2008), which 
requires robust theoretical models. 

Our results are exact (for the mass function) or are ex- 
pected to provide a good approximation (for the bias) in 
the limit of rare massive halos. However, this remains of 
interest as large-mass tails are also the most sensitive to 
cosmological parameters (e.g., through the linear growth 
factor and the primordial power spectrum), thanks to their 
steep dependence on mass or scale. Moreover, we think that 
more general fitting formulae (such as the one we provide 
for the mass function) should follow such theoretical pre- 
dictions in their relevant limits, so as to reduce the number 
of free parameters and improve their robustness. 
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